--- layout: page title: PO2 gradiants ---
import numpy as np
import pandas as pd
import statsmodels.api
import scipy.stats
import bokeh.io
import bokeh.plotting
#local .py file for some plotting functions and non-parametric bootstrapping utils
import plotting_utils
import numba
bokeh.io.output_notebook()
We can now read in the data into a dataframe for analyis.
df = pd.read_csv("./20190322_supp_table_2.csv")
We will add some calculations/generate a species averaged version of the dataframe.
df['species_underscore'] = [spec.replace(" ", "_") for spec in df['species']]
df_averages = df.groupby(['species', 'species_underscore', 'spiracle'], as_index=False).aggregate(np.average)
df_averages['subfamily'] = df.groupby(['species', 'species_underscore', 'spiracle'], as_index=False).aggregate(max)['subfamily']
species_per_subfam=df_averages.groupby(['subfamily', 'spiracle'], as_index=False).count().groupby('subfamily').aggregate(max).reset_index()[['subfamily', 'species']]
species_per_subfam.columns = ('subfamily', 'subfam_count')
df_averages = df_averages.merge(species_per_subfam, on='subfamily')
df_averages['log area (mm^2)'] = np.log10(df_averages['area (mm^2)'])
df_averages['log dist'] = np.log10(df_averages['depth (mm)'])
df_averages['log mass (g)'] = np.log10(df_averages['mass (g)'])
df_averages['log area/dist'] = np.log10(df_averages['area (mm^2)']/df_averages['depth (mm)'])
df_averages['log area^2/dist'] = np.log10(df_averages['area (mm^2)']**2/df_averages['depth (mm)'])
df['log area (mm^2)'] = np.log10(df['area (mm^2)'])
df['log dist'] = np.log10(df['depth (mm)'])
df['log mass (g)'] = np.log10(df['mass (g)'])
df['log area/dist'] = np.log10(df['area (mm^2)']/df['depth (mm)'])
df['log area^2/dist'] = np.log10(df['area (mm^2)']**2/df['depth (mm)'])
df_averages.head()
Using these morphological properties of the spiracles that we measured, we can calculate real world meaningful gas exchange properties with some unit conversions/physical gas properties. To calculate the diffusive capacity of the spiracles, we use the following relation:
\begin{align} G_{\mathrm{diff}} = \frac{\text{area}}{\text{depth}} * D * \beta \end{align}with $D$ (the diffusivity constant for O2 in air) $= 0.178 \text{ cm}^2$ and $\beta$ (the capacitance coefficient for oxygen in air) $= 404 \frac{\text{ nmol}}{\text{cm}^{3} \text{kPa}}$. With that value, you can calculate the oxygen consumption rate in $\frac{\text{nmol}}{\text{sec}}$ by multiplying the oxygen partial pressure gradient $\Delta \text{PO}_2$ across the spiracle in kPa by $G_{\text{diff}}$, that is
\begin{align} \text{Oxygen consumption rate }\left(\frac{\text{nmol}}{\text{sec}}\right) = \Delta \text{PO}_2 * G_{\mathrm{diff}} \end{align}To calculate advective gas transport capacities, we use $G_{\text{adv}}$, which comes from Poiseulle’s law and is
\begin{align} G_{\mathrm{adv}} = \frac{\text{area}^2}{8 * \pi * \text{dynamic viscosity} * \text{depth}} \end{align}where the dynamic viscosity of air is $1.86 * 10^{-8} \text{kPa sec}$. With these relations in hand, we can calculate both $G_{\text{diff}}$ and $G_{\text{adv}}$ for the beetle spiracles, both individually, and also by summing the individual capacities for all spiracles in the beetle.
df_test = df_averages.copy()
df_test['area/depth'] = (df_averages['area (mm^2)']/df_averages['depth (mm)'])/10
df_test['area^2/depth'] = ((df_averages['area (mm^2)']**2/df_averages['depth (mm)']))/(1000*1000*1000) #convert to cubic meters
df_test['g_diff'] = df_test['area/depth']*0.178*404
df_test['g_adv'] = df_test['area^2/depth']*(1/(np.pi*8*1.86*10**(-8)))
df_summed = df_test.groupby('species').median().reset_index()[['species', 'log mass (g)', 'g_diff', 'g_adv']]
df_summed['g_diff'] = df_test.groupby('species').sum().reset_index()['g_diff']
df_summed['g_adv'] = df_test.groupby('species').sum().reset_index()['g_adv']
plots = []
resample_size = 10_000
lw = 2
cs = 12
m = df_summed['log mass (g)']
g_diff = np.log10(df_summed['g_diff']*2) #double the number for g_diff since we only measured one side of the animal, they have two of each spiracle
slope, intercept = np.polyfit(m, g_diff, deg=1)
x = np.array([m.min(), m.max()])
y = slope * x + intercept
p = bokeh.plotting.figure(plot_height=300, plot_width=400,
x_range=(x[0]-0.1, x[1]+0.1), y_range=(g_diff.min()-0.1, g_diff.max()+0.2))
p.outline_line_color = None
p.yaxis.minor_tick_line_color = None
p.xaxis.minor_tick_line_color = None
p.grid.grid_line_color = None
slope_comp = 0.75
intercept1 = plotting_utils.first_intercept(slope_comp, x.max(), g_diff.min()) -0.5
line_scale = (y.max() - y.min())/5
around_line=0.2
for i in line_scale*np.array(range(30))+intercept1:
try:
lx, ly = plotting_utils.generate_line(intercept=i, slope=slope_comp,
bounds=(x[0]-around_line, x[1]+around_line,
g_diff.min()-around_line, g_diff.max()+around_line+0.4), point=x[1])
p.line(lx, ly, color='grey', alpha=0.3)
except:
pass
bs_slope_reps, bs_intercept_reps, _, _ = plotting_utils.draw_bs_pairs_linreg(m, g_diff, size=resample_size)
p.title.text = ('Slope: '+ str(round(slope, 2)) + ' intercept: ' + str(round(intercept, 2)) +' slope 95% CI: ' +
str([round(j, 3) for j in np.percentile(bs_slope_reps, [2.5, 97.5])]) + ' ' + str(np.sum(bs_slope_reps > 0.75)))
x_boot = np.linspace(m.min(), m.max(), 200)
y_boot = np.outer(bs_slope_reps, x_boot) + np.stack([bs_intercept_reps]*200, axis=1)
low, high = np.percentile(y_boot, [2.5, 97.5], axis=0)
p1 = np.append(x_boot, x_boot[::-1])
p2 = np.append(low, high[::-1])
p.patch(p1, p2, alpha=0.5, color='lightgrey')
p.circle(m, g_diff, color='black', size=cs)
p.line(x, y, color='black', line_width=lw, line_cap='round')
p.output_backend='svg'
plots.append(p)
#bokeh.io.show(p)
m = df_summed['log mass (g)']
g_adv = np.log10(df_summed['g_adv']*2) #double the number for g_diff since we only measured one side of the animal, they have two of each spiracle
slope, intercept = np.polyfit(m, g_adv, deg=1)
x = np.array([m.min(), m.max()])
y = slope * x + intercept
p = bokeh.plotting.figure(plot_height=300, plot_width=400, x_range=(x[0]-0.1, x[1]+0.1), y_range=(g_adv.min()-0.2, g_adv.max()+0.5))
p.outline_line_color = None
p.yaxis.minor_tick_line_color = None
p.xaxis.minor_tick_line_color = None
p.grid.grid_line_color = None
slope_comp = 0.75
intercept1 = plotting_utils.first_intercept(slope_comp, x.max(), g_adv.min())
line_scale = (y.max() - y.min())/7
around_line=0.4
for i in line_scale*np.array(range(30))+intercept1:
try:
lx, ly = plotting_utils.generate_line(intercept=i, slope=slope_comp,
bounds=(x[0]-around_line, x[1]+around_line,
g_adv.min()-around_line, g_adv.max()+around_line+0.5), point=x[1])
p.line(lx, ly, color='grey', alpha=0.3)
except:
pass
bs_slope_reps, bs_intercept_reps, _, _ = plotting_utils.draw_bs_pairs_linreg(m, g_adv, size=resample_size)
p.title.text = ('Slope: '+ str(round(slope, 2)) + ' intercept: ' + str(round(intercept, 2)) +' slope 95% CI: ' +
str([round(j, 3) for j in np.percentile(bs_slope_reps, [2.5, 97.5])]) + ' ' + str(np.sum(bs_slope_reps < 0.75)))
x_boot = np.linspace(m.min(), m.max(), 200)
y_boot = np.outer(bs_slope_reps, x_boot) + np.stack([bs_intercept_reps]*200, axis=1)
low, high = np.percentile(y_boot, [2.5, 97.5], axis=0)
p1 = np.append(x_boot, x_boot[::-1])
p2 = np.append(low, high[::-1])
p.patch(p1, p2, alpha=0.5, color='lightgrey')
p.circle(m, g_adv, color='black', size=cs)
p.line(x, y, color='black', line_width=lw, line_cap='round')
p.output_backend='svg'
plots.append(p)
bokeh.io.show(plots[0])
bokeh.io.show(plots[1])
To calculate the ∆PO2 across the spiracles needed to supply beetle metabolic demand by diffusion, the metabolic rate for a quiescent beetle at a body temperature of 25 °C of a given mass was estimated from here with the following equation: $\text{log}_{10} (\text{metabolic rate } \mu \text{W}) = 3.20 + 0.75*\mathrm{log}_{10}(\text{mass (g)}$ and the assumption of 20.7 kJ/L for O2. For O2 at 25 °C, 24.5 mol/L was also assumed. Based on here, flight metabolic rate of small insects is in the range of 8x resting, whereas it is on the order of 32x resting metabolic rate in large insects. Most scarab beetles are endothermic during flight, so flight metabolic rates of these warm beetles could double this value to 64x with a 10 °C increase in thoracic temperature if the Q10 is 2. Even this calculation may be a conservative estimate of hovering metabolic rate, as oxygen consumption rose to 90x higher than those of quiescent, 25 °C fig beetles in one study, and maximal flight metabolic rates may be 1.25-2x higher than during hovering. Therefore, we estimated maximal aerobic metabolic rate during flight as 90x those of quiescent beetles. The required PO2 gradient across the spiracles to support gas exchange by diffusion at rest and during flight was calculated by rearranging our previous $G_{text{diff}}$ equation and performing unit conversions as follows:
With this relation in hand, we plot the required oxygen partial pressure gradient across just the spiracle needed to supply the metabolic demand of the insect.
m = df_summed['log mass (g)']
g_diff = np.log10(df_summed['g_diff']*2)
slope, intercept = np.polyfit(m, g_diff, deg=1)
po2 = np.log10((1*((10**(3.20 + 0.75*m))*(1/20.7)*(1/24.5))/(df_summed['g_diff']*2)))
print(slope)
N = 200
y_top = 90
x = np.linspace(m.min(), m.max(), N)
y = np.linspace(1, y_top, N)
im = np.zeros((N, N))
for j, yi in zip(range(N), y):
for i, xi in zip(range(N), x):
im[i, j] = (yi*((10**(3.20 + 0.75*xi))*(1/20.7)*(1/24.5))/(10**(slope * xi + intercept)))
N = 800
y_top = 90
x = np.linspace(m.min(), m.max(), N)
y = np.linspace(1, y_top, N)
x_pos = []
for j, yi in zip(range(N), y):
[x_pos.append([xj, yj, zj]) for xj, yj, zj in zip(x, np.log10(yi*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/(10**(slope * x + intercept))),
yi*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/(10**(slope * x + intercept)))]
x_min, y_min, v_min = (np.array(x_pos)[:, 0].min(), np.array(x_pos)[:, 1].min(), np.array(x_pos)[:, 2].min())
x_stride, y_stride = ((np.array(x_pos)[:, 0].max() - x_min)/N, (np.array(x_pos)[:, 1].max() - y_min)/N)
im = np.ones((N, N))*v_min
for xj, yj, vj in x_pos:
im[int(np.ceil((yj-y_min)/y_stride-1)), int(np.ceil((xj-x_min)/x_stride-1))] = vj
im = scipy.ndimage.gaussian_filter(im, sigma=1)
p = bokeh.plotting.figure(tooltips=[("x", "$x"), ("y", "$y"), ("value", "@image")], plot_height=300, plot_width=400)
#p.x_range.range_padding = p.y_range.range_padding = 0
cmap = bokeh.models.LinearColorMapper(palette='Viridis256', low=im.min(), high=im.max())
cmap_low, cmap_high = (np.min(1*(10**po2)), np.max(90*(10**po2)))
cmap = bokeh.models.LinearColorMapper(palette='Viridis256', low=cmap_low, high=cmap_high)
p.image(image=[im], x=x.min(), y=np.log10(im.min()), dw=x.max()-x.min(), dh=np.log10(im.max())-np.log10(im.min()), color_mapper=cmap, level="image", )
color_bar = bokeh.models.ColorBar(color_mapper=cmap, location=(0,0), ticker=bokeh.models.BasicTicker(desired_num_ticks=12, base=10))
p.add_layout(color_bar, 'right')
p.grid.grid_line_color = None
p.patch([x.min(), x.max(), x.max(), x.min()], [np.log10((8*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).min(),
np.log10((8*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).max(),
np.log10(im.min()), np.log10(im.min())], color='white', line_width=2)
p.patch([x.min(), x.max(), x.max(), x.min()], [np.log10((90*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).min(),
np.log10((90*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).max(),
np.log10(im.max()), np.log10(im.max())], color='white', line_width=2)
p.line([x.min(), x.max()], [np.log10(21), np.log10(21)], color='lightgrey', line_width=2, alpha=0.75)
#p.patch([x.min(), x.max(), x.max(), x.min()], [np.log10((1*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).min(),
# np.log10((1*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).max(),
# np.log10((90*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).max(),
# np.log10((90*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).min()], color='white', line_width=0, alpha=0.2)
line_color='black'
#[p.line(x, np.log10((yi*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))), color=line_color) for yi in np.linspace(1, y_top, 10)]
[p.line(x, np.log10((yi*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))), color=line_color, line_width=lw, line_cap='round') for yi in [1, 8, 90]]
#[p.line(x, np.log10((yi*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(0.33 * x + intercept))), color='white', line_dash='dashed', line_width=5, line_cap='round') for yi in [1, 8, 90]]
dot_color='white'
dot_line='black'
dot_width = 1
cmapper = bokeh.transform.linear_cmap('c', palette='Viridis256', low=np.min(1*(10**po2)), high=np.max(90*(10**po2)))
for mult in [1, 8, 90]:
source = bokeh.models.ColumnDataSource(data=dict(x=m, y=np.log10(mult*(10**po2)), c=mult*(10**po2)))
p.circle('x', 'y', source = source, size=cs, line_color=dot_line, line_width = dot_width, fill_color=cmapper)
p.outline_line_color = None
p.yaxis.minor_tick_line_color = None
p.xaxis.minor_tick_line_color = None
p.output_backend='svg'
plots.append(p)
bokeh.io.show(bokeh.layouts.gridplot([plots[0], plots[2], plots[1]], ncols=3))
The required partial pressure gradient with a 90x resting metabolic demand is
90*(10**po2)
These values exceed the partial pressure of air in the atmosphere! No way a large beetle can use diffusion for gas exchange, it must use active advective processes.
#TO EXPORT THE PLOTS FROM ABOVE
#bokeh.io.export_svgs(plots[2], filename="./ptest.svg")
%reload_ext watermark
%watermark -p bokeh